#%% [markdown]
"""
# 时间需预测问题
 
### 作业流程
+ 先导入数据,并对数据进行转换
+ 绘制时间序列图(折线图)
+ 观察3年销售额时间序图,总结
+ 预测第四年 ~~(ARIMA)~~ (winter模型) 
+ 使用的模块
    + pandas 表格处理
    + matplotlib 基础绘图
    + statsmodels 统计模块(自相关图,偏相关图等)
    + scipy.stats

扩展阅读
https://cloud.tencent.com/developer/article/1646121
https://blog.csdn.net/sanjianjixiang/article/details/103939699
"""
#%%
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf  # 非常规库,如果没有请剔除
from statsmodels.tsa.stattools import adfuller

# from statsmodels.tsa.arima.model import ARIMA  # ARIMA模型导入
from scipy import stats

plt.rcParams["font.family"] = "Microsoft YaHei"
# %% [markdown]
""" 
### 0.数据预处理
+ 导入数据
+ 行列转换,将数据从12行3列转为36行2列
"""
#%%
df = pd.read_excel("作业数据.xlsx", sheet_name=1)
df.columns = df.iloc[0, :]
df = df.iloc[1:, :]
# 对行列进行转换
con = []
for y_name in df.columns[1:]:
    tdf = df.loc[:, ["月份", y_name]].rename(columns={y_name: "销售额"})
    tdf["年"] = y_name
    tdf["年-月"] = tdf["年"] + "-" + tdf["月份"].astype(str)
    tdf["销售额"] = tdf["销售额"].astype(float)
    con.append(tdf)
tdf = pd.concat(con)
tdf = tdf.reset_index().drop(columns="index")
tdf.info()
display(tdf)
#%% [markdown]
""" 
### 1.展示过去三年销售额时间序列图
+ 折线图绘制
绘制后发现,数据中存在较强的周其性,但是趋势不显著
+ 趋势获取,
先不对数据进行处理
"""
#%%
# 直接计算趋势线
linear = np.polyfit(tdf.index, tdf["销售额"], deg=1)
y_prodict = np.polyval(linear, tdf.index)

#%% [markdown]
"""  
依照模型发现,整体趋势成持平可以看到模型中截距为 $a_1=-0.01$
- 但是由于没有剔除异常波动,需要对其进行周期计算
"""
#%%
plt.figure(figsize=(12, 5))
plt.plot(tdf.index, tdf["销售额"], label="销售额", marker="o")
plt.plot(tdf.index, y_prodict, label="趋势线", linestyle=":", c="r")
plt.xticks(tdf.index, tdf["年-月"].values, rotation=-90)
plt.text(
    tdf.index[0],
    y_prodict[0],
    s="直接计算趋势y={:.2f}x+{:.2f}".format(linear[0], linear[1]),
    color="r",
)

#%% [markdown]
"""
### 2.观察过去三年销售额时间序列图，总结销售额变化特点和规律
- 由之前模型总结得出
    - 1. 数据存在一定的周期性,目测为12个月一周期
    - 2. 整体趋势成轻微下滑,但是由于存在一定的扰动,需要剔除
- 结论
    - 通过ACF得知,数据以12个单位为一个周期,并且成拖尾
    - 通过PACF得知,模型趋势性较弱
    - ADF检验通过,得知模型相对平稳无需差分
- ARIMA选择AR模型
"""
# %% [markdown]
"""
#### 自相关系数(ACF)
- 不难看出,自相关系数(ACF)中呈现了拖尾,且具有较好的周期性
- 从第一个节点0点处向后,因存在2个余弦周期,(cos(pi*x)/2) 第一个落点为3
- 通过图像观察,第一个完整周期也出现在12附近
"""

"""
#### 偏自相关系数(PACF) 严格自相关系数
- 呈现截尾
"""
#%%
# 另一个比较粗暴的解决方法
# _tdf = tdf.copy()
# for i in range(30):
#     _tdf["销售额:" + tdf["年-月"].iloc[i]] = _tdf["销售额"].shift(i + 1).astype(float)

# corrs = _tdf.corr()
# plt.figure(figsize=(12, 12))
# plt.imshow(corrs, cmap=plt.cm.RdBu_r)

# self_corr = _tdf.corr()["销售额"]
# plt.figure(figsize=(12, 5))
# plt.bar(self_corr.index, self_corr.values)
# plt.xticks(self_corr.index, rotation=-90)
#%%
# 1. 周期检查 acf 与 pacf
fig = plt.figure(figsize=(12, 10))
ax1 = fig.add_subplot(2, 1, 1)
plot_acf(tdf["销售额"].values, ax=ax1)

ax2 = fig.add_subplot(2, 1, 2)
plot_pacf(tdf["销售额"].values, ax=ax2)

# %% [markdown]
"""
#### CMA 和 差分检验
为了验证周期情况通过,通过CMA(中心移动平均发进行检验)和ADF进行检验
adf_api: https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.adfuller.html
+ adf 检验后,adf值
"""
#%%
adf_res = adfuller(
    tdf["销售额"],
    maxlag=13,
    regression="ct",  # 常数+趋势
)
display(
    dict(
        zip(["adf值", "置信度p", "所使用的之后数", "观测值数量", "临界值", "icbest", "resstore"], adf_res)
    )
)

# %%
# CMA均值也与adf检验结果一致,发现已12为周期波动较小
CMA = tdf["销售额"].rolling(12).mean().shift(-6).dropna()

#%%
linear = np.polyfit(CMA.index, CMA.values, deg=1)
cma_prodict = np.polyval(linear, tdf.index)

plt.figure(figsize=(12, 5))
plt.plot(tdf.index, tdf["销售额"], label="销售额", marker="o")
plt.plot(tdf.index, cma_prodict, label="CMA", linestyle=":", c="r")
plt.xticks(tdf.index, tdf["年-月"].values, rotation=-90)


plt.text(
    tdf.index[0],
    cma_prodict[0],
    s="直接计算趋势y={:.2f}x+{:.2f}".format(linear[0], linear[1]),
    color="r",
)
#%%
# 验证模型准确度
X = sm.add_constant(CMA.index.tolist())
ols_model = sm.OLS(CMA.values.tolist(), X, hasconst=False).fit()
display(ols_model.summary())
#%% [markdown]
"""
+ 可以发现,模型拟合效果极佳,$R^2$达到了1
+ F检验和T检验p值显著性均通过检验(在正态分布的情况下,Y'和Y方差一致,x值均值一致)两者均通过显著性检验,模型有效
"""
#%% [markdown]
""" 
### 3.预测第四年1-12月的销售额
通过之前检验
"""
#%%
# ARIMA
# ar = ARIMA(tdf["销售额"],order=(10,0,0)).fit()
# print(ar.summary())
# plt.plot(
# ar.forecast(12, alpha=0.05))

# %% [markdown]
""" 
+ 计算季节因子,通过 $ 销售额 / CMA $
"""
#%%
tdf = pd.merge(
    tdf,
    CMA,
    how="left",
    left_index=True,
    right_index=True,
)
tdf["secson_index"] = tdf["销售额_x"] / tdf["销售额_y"]
#%%
secson_index = tdf[["月份", "secson_index"]].groupby(by=["月份"]).mean()  # 计算各季节指数
display(secson_index)
display("季节指数均值{}".format(secson_index.secson_index.mean()))
#%%
# 季节指数校准
secson_index = secson_index / secson_index.secson_index.mean()
# 重新赋值
tdf = pd.merge(tdf, secson_index, how="left", left_on="月份", right_index=True)
# 进行预测
t = []
for n in range(4):
    t.extend([i + 1 for i in range(12)])
t
_df = pd.DataFrame(t, columns=["month"])
_df = pd.merge(_df, secson_index, how="left", left_on="month", right_index=True)
_df["sells"] = (_df.index * ols_model.params[1] + ols_model.params[0]) * _df[
    "secson_index"
]

_df["year"] = _df.index // 12 + 1
_df["year-month"] = "y:" + _df["year"].astype(str) + "-m:" + _df["month"].astype(str)
_df
#%% [markdown]
""" 
#### 结果假设下检验
+ 配对t检验: 双尾$ p=0.532552969007966 $ 远大于 $\\a = 0.05 $
    + 接受 $ H_0 $ 预测值与实际值之间平均值不存在显著差异
+ ANOVA(方差一致性检验):$ p=0.9758031118109965 $ 远大于 $ \\a = 0.05 $
    + $ H_0 $ 预测值与实际值不存在显著差异,二类错误,故接受原假设
"""
#%%
display(stats.ttest_rel(tdf["销售额_x"], _df.loc[: tdf.shape[0] - 1, "sells"]))
display(stats.f_oneway(tdf["销售额_x"], _df.loc[: tdf.shape[0] - 1, "sells"]))
display(
    "两组样本的平均误差为{}".format(
        (sum((tdf["销售额_x"] - _df.loc[: tdf.shape[0] - 1, "sells"]) ** 2) / tdf.shape[0])
        ** (1 / 2)
    )
)

# %% [markdown]
""" 
+ 通过模型得出,
    + 整体收入依然成缓慢增长,每年增长率约为 $1.07^{12}$
    + 初始截距为169.94
"""
#%%
plt.figure(figsize=(12, 5))
plt.plot(tdf.index, tdf["销售额_x"], label="原销售额", marker="o")
plt.plot(_df.index, _df.sells, label="预测销售额", linestyle=":", c="r")
plt.xticks(_df.index, _df["year-month"].values, rotation=-90)
plt.plot([36], [295], label="要求参考点", marker="o", c="g")
plt.text(
    _df.index.tolist()[0],
    _df.sells[0] + 30,
    s="直接计算趋势y={:.2f}*x*secson_index+{:.2f}".format(
        ols_model.params[1], ols_model.params[0]
    ),
    color="r",
)
plt.legend()
# %% [markdown]
""" 
### 假设第四年1月份销售额为295万元，你的预测误差为多少？
"""
# %%
print("第四年第一个月,预测值与实际值的误差为{:.5}万元".format(abs(_df.loc[36, "sells"] - 295)))
# %%
